# -*- coding: utf-8 -*-
"""
Created on Mon Aug 23 17:46:16 2021

@author: admin
"""


import os
from osgeo import gdal
import sys



def ClipDEMByArea(OutFilePath,DEMFilePath,SelectArea):
    '''
    根据选择范围裁剪DEM,并输出
    agrs:
        outFilePath:裁剪DEM输出地址
        DEMFilePath:被裁减DEM地址
        SelectArea:list  [(xmin,ymax),(xmax,ymin)] 框选范围 左上角，右下角
    '''
    DEM_ptr=gdal.Open(DEMFilePath) 
    DEM_GeoTransform=DEM_ptr.GetGeoTransform() # 读取影像的投影变换
    DEM_InvGeoTransform=gdal.InvGeoTransform(DEM_GeoTransform)
    SelectAreaArrayPoints=[gdal.ApplyGeoTransform(DEM_InvGeoTransform,p[0],p[1])  for p in SelectArea]
    SelectAreaArrayPoints=list(map(lambda p:(int(p[0]),int(p[1])),SelectAreaArrayPoints)) # 确定坐标
    
    [(ulx,uly),(brx,bry)]=SelectAreaArrayPoints
    rowCount,colCount=bry-uly,brx-ulx
    
    # 输出DEM的桌面坐标转换
    Out_Transfrom=list(DEM_GeoTransform)
    Out_Transfrom[0]=SelectArea[0][0]
    Out_Transfrom[3]=SelectArea[0][1]
    
    # 构建输出DEM
    Bands_num=DEM_ptr.RasterCount
    gtiff_driver = gdal.GetDriverByName('GTiff')
    out_dem = gtiff_driver.Create(OutFilePath, colCount, rowCount, Bands_num)
    out_dem.SetProjection(DEM_ptr.GetProjection())
    out_dem.SetGeoTransform(Out_Transfrom)
    
    for i in range(1,Bands_num+1):
        data_band=DEM_ptr.GetRasterBand(i)
        out_band=out_dem.GetRasterBand(i)
        data=data_band.ReadAsArray(ulx, uly, colCount, rowCount)
        out_band.WriteArray(data)
    del out_dem
        

if __name__=="__main__":
    outFilePath=r"/mnt/d/Orthorectification/data/clip2.tif"
    DEMFilePath=r'/mnt/d/Orthorectification/data/mosaic9(1).tif'
    Selectarea=[(112.5,24),(113.5,23)]
    ClipDEMByArea(outFilePath,DEMFilePath,Selectarea)
    pass


